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SECTION  1 
INTRODUCTION 


The  prediction  of  blast  wave  environments  via  numerical  simulation  is  one  of  the 
organizing  themes  of  the  work  discussed  here  and  it’s  original  motivation.  Such 
environments  can  be  categorized  as  ideal  or  nonideal,  as  well  as  early  time  (equiv¬ 
alently,  high  overpressure)  or  late  time  (equivalently,  low  overpressure).  The  work 
reported  on  here  is  largely  based  on  code  technology  developed  in  our  preceding 
efforts  at  modeling  high  overpressure  environments,  both  ideal  and  nonideal  (using 
thermal  layer  models)  [17].  Our  code  development  efforts  over  the  past  few  years 
have  been  directed  at  extending  this  technology  in  the  following  directions:  (1) 
problems  involving  wave  number  stiffness  which  include  as  a  special  case  nearly  in¬ 
compressible  flow,  (2)  problems  for  which  the  incompressible  equations  themselves 
is  appropriate,  (3)  problems  involving  source  term  stiffness  and  (4)  problems  in¬ 
volving  strong  shock  wave  propagation  in  materials  with  complicated  constitutive 
properties.  The  applications  which  have  focused  this  work  arise  largely  from  very 
late  time  blast  wave  environments,  especially  the  problems  inherent  in  the  calcula¬ 
tion  of  fireball  rise.  The  motivation  for  our  work  in  area  (4)  above  comes,  instead, 
from  a  desire  to  understand  the  effects  of  high  overpressure  blast  waves  on,  for 
example,  hardened  targets  or  differing  soil  types. 

If  one  considers  the  problem  of  fireball  rise,  it  is  immediately  apparent  that  the 
analogue  computational  problem  is  really  at  least  three  problems:  it’s  initial  for¬ 
mation,  a  long  intermediate  stage  in  which  it  begins  to  rise  from  a  nearly  stationary 
position  and  a  late  stage  in  which  the  local  Mach  number  gradually  increases  and 
compressibility  effects  are  once  again  important.  For  our  purposes  here,  we  con¬ 
sider  the  first  part  as  solved  by  our  older  work.  The  second  stage  can  be  solved 
either  by  considering  an  incompressible  model  and  solving  the  derived  equations  or 
by  developing  specialized  implicit  strategies  in  the  context  of  a  compressible  flow 
solver.  There  is  a  stronger  requirement  for  the  latter  approach  as  the  Mach  number 
increases  at  later  times.  On  the  other  hand,  related  problems  such  as  the  prediction 
of  contaminant  transport  in  a  battlefield  or  theatre  environment  axe  better  handled 
by  using  an  appropriate  model  with  sonic  waves  eliminated  a  priori.  Accordingly, 
we  have  worked  on  both  approaches.  At  the  same  time,  we  have  especially  concen¬ 
trated  on  particulate  transport  and  reactive  flow  in  these  flow  fields,  since  they  axe 
critically  important  inputs  (and  outputs)  for  the  intended  applications.  This  has 
led  to  our  work  on  source  term  stiffness. 
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1.1  IMPLICIT-EXPLICIT  SCHEMES. 


In  our  report  [17],  we  discussed  the  early  stages  of  the  development  of  the  Implicit- 
Explicit  (I-E)  scheme,  in  particular  some  one-dimensional  inviscid  and  viscous  flow 
calculations,  see  [10]— [11],  and  we  also  presented  future  tasks.  Several  studies  have 
been  completed  since  and  they  represent  major  advances  in  this  area.  We  start 
with  a  brief  description  of  these  studies.  First,  Collins  has  successfully  accom¬ 
plished  several  objectives  that  extended  the  generic  I-E  methodology;  these  were 
documented  in  great  detail  in  his  Ph.D  dissertation  [12].  Some  aspects  of  this  work 
were  presented  at  the  First  European  Computational  Fluid  Dynamics  conference, 
see  [14].  A  summary  of  the  objectives  and  results  of  the  work  carried  out  by  Collins 
[12]  and  Collins  et.  al.  [14]  are  given  in  Section  1.1.1  below.  Second,  the  I-E 
methodology  was  implemented  for  flows  that  are  modeled  by  systems  with  source 
terms,  in  particular  chemically  reacting  and  dusty  gas  models,  including  the  cases 
of  stiff  sources,  see  [19]  and  [22].  The  results  of  these  studies  are  summarized  in 
Section  1.1.2.  below.  Third,  an  early  version  of  the  work  by  Collins  et.  al  (in  the 
context  of  ID  flows)  was  submitted  to  and  accepted  for  publication  recently  [13]. 

1.1.1  Zero  Mach  Number  Limit. 

The  main  objective  of  the  work  done  by  Collins  [12]  and  Collins  et  al  [13]  [14]  is  to 
eliminate  the  explicit  Courant-Friedrichs-Levy  (CFL)  constraint  on  unsteady  com¬ 
pressible  flow  calculations  by  using  a  high  resolution  explicit  scheme  in  conjunction 
with  an  implicit  scheme  in  such  a  way  that  the  explicit  CFL  restriction  is  miti¬ 
gated  and  that  the  solution  retains  much  of  the  temporal  accuracy  provided  by  the 
explicit  scheme.  Explosion  dynamics  provides  a  particularly  challenging  example 
where  this  methodology  is  needed.  There  are  two  phases  of  interest  in  these  types  of 
problems:  the  initial  blast  wave  and  the  subsequent  cloud  rise.  Typically,  the  blast 
wave  phase  is  computed  using  a  high  resolution  explicit  compressible  code  while  the 
cloud  rise  is  simulated  with  an  incompressible  code.  At  late  time  during  the  cloud 
rise  phase,  the  Mach  number  becomes  significant  once  again  and  can  approach  one; 
thus  violating  the  assumption  of  incompressibility.  The  I-E  strategy  could  be  used 
to  handle  both  phases  of  the  problem,  and  in  addition,  it  offers  significant  savings 
for  just  the  blast  phase  of  this  problem.  For  example,  if  the  temperature  in  the 
fireball  is  too  high  and  the  explicit  time  step  of  such  calculation  becomes  too  small 
for  such  calculation  to  be  feasible,  the  scheme  would  allow  for  a  much  larger  time 
step  without  significant  degradation  in  the  solution.  At  late  times,  the  same  code 
can  be  used  to  compute  the  cloud  rise.  As  the  Mach  number  in  the  cloud  increases, 
the  I-E  method  will  transition  continuously  to  a  high  resolution  explicit  scheme. 

The  study  followed  the  guidelines  stated  below: 

(1)  Develop  an  Eulerian  implicit-explicit  scheme,  similar  to  the  one  in  reference 
[16],  which  satisfies  several  design  principles  (for  details,  see,  e.g.,  [14]). 

(2)  Use  the  second-order  Godunov  scheme  as  the  underlying  high  resolution 
explicit  scheme. 
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(3)  Extend  this  method  to  two  space  dimensions  using  the  unsplit  Godunov 
scheme  in  [7]. 

(4)  Provide  numerical  calculations  which  demonstrate  the  feasibility  of  this 
method. 

The  work  consists  of  two  parts.  In  Part  I,  the  method  is  implemented  for  the 
special  case  of  one-dimensional  inviscid  compressible  flow  where,  for  convenience, 
a  polytropic  equation-of-state  is  used.  The  methodology  can  be  easily  extended 
to  gas  dynamics  with  a  general  equation-of-state  using  the  ideas  in  [8].  A  more 
difficult  question  is  the  extension  of  these  ideas  to  two  or  three  space  dimensions; 
this  was  done  in  Part  II  for  the  two-dimensional  case.  Implementations  for  systems 
in  higher  space  dimensions  follows  directly  from  the  two-dimensional  case. 

At  first,  the  implicit-explicit  ideas  are  developed  for  the  special  case  of  linear  ad- 
vection  in  one  space  dimension;  extensions  of  this  linear  scalar  scheme  to  handle 
source  terms  follows.  In  particular,  it  is  shown  that  this  scheme  satisfies  the  design 
principles,  thereby  fixing  these  ideas  for  extension  to  nonlinear  systems.  Also,  the 
potential  for  nonlinear  instabilities  is  recognized  and  an  indication  is  given  as  to 
how  to  handle  such  instabilities.  A  nonlinear  stability  constraint  is  introduced  for 
the  method  in  the  special  case  of  a  scalar  conservation  law.  Next,  the  scheme  is 
described  for  systems  of  conservation  laws;  the  formulae  necessary  for  compressible 
gas  dynamics  are  included.  For  that  case,  it  is  assumed  that  a  sufficiently  smooth 
numerical  flux  function  exists.  Since  Newton’s  method  is  used  to  solve  the  system 
in  one  space  dimension,  a  key  component  of  the  work  in  Part  I  is  the  introduction 
of  an  appropriately  smooth  numerical  flux  function  for  the  hyperbolic  equations.  A 
suitable  version  of  the  Engquist-Osher  flux  [6]  is  used,  as  modified  for  systems  by 
Bell  et.  al.  [4].  The  version  used  is  sufficiently  smooth  so  that  the  Newton’s  method 
linearization  is  well-behaved;  in  particular,  it  converges  to  steady  states  even  in  the 
presence  of  strong  shocks.  The  construction  of  this  flux  function,  together  with 
several  approximations  that  are  efficient  for  practical  large-scale  applications  are 
presented  in  detail.  For  the  two-dimensional  system,  an  approximation  to  the 
Godunov  flux  is  used  and  a  nonlinear  multigrid  technique  is  applied  to  solve  the 
resulting  equations.  To  conclude  the  first  part,  numerical  results  are  presented 
for  the  one-dimensional  case.  It  is  shown  that  the  schemes  are  able  to  propagate 
slowly  moving  shocks  and  contact  surfaces  through  a  Cartesian  mesh.  Also,  the 
scheme’s  ability  to  capture  a  stationary  shock  in  a  diverging  duct  is  demonstrated. 
A  numerical  study  is  presented  which  demonstrates  second-order  accurate  steady 
states. 

In  Part  II,  the  scheme  is  generalized  to  two  dimensions  using  Colella’s  Corner 
Transport  Upwind  (CTU)  scheme  [7]  as  the  underlying  explicit  scheme.  Once 
again,  the  focus  is  on  the  scalar  linear  advection  equation,  but  this  time  in  two  space 
dimensions.  The  stability  results  here  are  far  different  from  the  one-dimensional 
case.  The  scheme  of  choice,  i.e.,  the  one  studied  in  Part  I,  is  unstable  when  applied 
to  the  scalar  linear  advection  in  two  dimensions;  therefore,  in  an  effort  to  find  a 


3 


better  implicit  scheme,  a  two  parameter  family  of  schemes  is  considered,  with  the 
same  simple  structure  as  the  implicit  scheme  in  Part  I.  A  stable  scheme  is  found 
but  continuity  in  CFL  number  is  lost.  It  is  shown  that  even  the  stable  implicit 
scheme  is  unstable  when  applied  to  the  scalar  linear  advection  equation  in  a  hybrid 
mode.  Because  of  this,  the  focus  of  attention  shifts  to  the  nearly  incompressible 
flow  application  under  the  constraints  that  the  particle  characteristics  are  always 
handled  explicitly  (that  is  the  advection  part  of  the  flow)  whereas  only  the  acoustic 
modes  axe  treated  implicitly. 

Next,  the  truncation  error  and  stability  properties  of  the  hybrid  schemes  are  investi¬ 
gated,  as  applied  to  a  linear  system  derived  from  the  isentropic  gas  dynamic  condi¬ 
tion  equations.  The  analysis  is  done  under  the  assumption  of  nearly  incompressible 
flow.  The  stability  analysis  indicates  that,  under  not  very  restrictive  constraints  on 
the  advection  part  of  the  flow,  the  original  scheme  is  stable.  Very  similar  results 
are  obtained  with  the  unconditionally  stable  implicit  scheme  discussed  earlier.  The 
truncation  error  analysis  verifies  the  order  of  accuracy  for  the  hybrid  scheme,  and 
in  particular,  verifies  second  order  accurate  steady  states.  It  also  shows  why  one 
cannot  expect  the  scheme  to  work  at  arbitrarily  low  Mach  numbers,  problem  for 
all  compressible  codes  [29]. 

Part  II  is  concluded  with  a  presentation  of  computational  results  for  two  problems: 
a  doubly  periodic  shear  layer  and  a  simulation  of  one  of  the  Brown  and  Roshko 
shear  layer  experiments  [5].  Both  of  these  are  examples  of  unstable,  nearly  incom¬ 
pressible  flow.  These  results  show  that  the  scheme  performs  as  desired;  however, 
the  efficiency  of  the  nonlinear  multigrid  algorithm  used  in  this  study  (and  described 
in  detail  in  [14])  leaves  room  for  improvement. 

To  summarize  this  work  we  examine  the  results.  The  results  of  the  ID  cases  show 
that  explicit  high  resolution  upwind  schemes  can  be  smoothly  hybridized  on  a 
mode-by-mode  basis  with  an  implicit  method  so  that  problems  involving  localized 
wave  speed  stiffness  and  convergence  to  steady  state  can  be  solved  without  sacri¬ 
ficing  accuracy  and  resolution  for  strong  wave  interactions  and  energy  containing 
modes.  Several  of  the  calculations  illustrate  the  importance  of  replacing  the  Go¬ 
dunov  flux  by  the  smoother  approximate  Engquist-Osher  flux  when  the  system 
is  linearized.  Questions  involving  the  optimal  choice  of  an  approximate  Riemann 
solver  and  appropriate  dissipative  mechanisms  have  not  been  completely  resolved; 
however,  a  very  simple  and  efficient  solver  augmented  by  a  small  amount  of  artifi¬ 
cial  viscosity  and  some  additional  dissipation  at  sonic  points  has  been  capable  of 
handling  our  test  problems. 

The  stability  results  for  this  methodology  in  two  space  dimensions  forced  us  to 
focus  Part  II  on  the  low  speed  nearly  incompressible  flow  application.  The  com¬ 
putational  results  for  this  part  demonstrate  that  the  implicit-explicit  methodology 
is  well  suited  for  this  application.  The  results  show  that  for  low  Mach  number 
flows,  the  degradation  of  the  solution  due  to  the  implicit  part  of  the  computation 
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is  minimal  and  the  high  resolution  capability  of  the  underlying  explicit  scheme  is 
preserved.  Also,  since  both  the  implicit  and  explicit  results  compare  well  with 
the  incompressible  results  reported  in  reference  [3],  the  approximate  phase  space 
solution  and  the  Godunov  flux  discussed  in  detail  in  [14]  appear  to  be  adequate. 
The  major  drawback  with  the  present  implementation  of  this  scheme  is  the  cost  of 
solving  the  nonlinear  system.  For  the  current  multigrid  implementation,  calcula¬ 
tions  must  be  executed  at  a  Courant  number  of  approximately  7  or  higher  to  be 
competitive  with  the  explicit  scheme;  this  corresponds  to  Mach  numbers  of  0.1  or 
lower.  However,  Moo  <  -1,  it  is  questionable  which  solution  strategy  is  competitive 
with  the  explicit  code. 

We  remark  that  an  attempt  was  made  to  accelerate  the  convergence  of  the  resulting 
nonlinear  system  of  equations  by  following  the  preconditioning  matrix  idea  of  Shuen 
et  al  [27].  This  idea  consists  of  introducing  an  extra  derivative  with  respect  to 
some  pseudo-time;  this  additional  finite  difference  derivative  is  premultiplied  by 
preconditioning  matrix  and  in  the  limit  of  steady-state  with  respect  to  pseudo¬ 
time  the  solution  to  the  new  system  converges  to  the  solution  of  the  original  (real) 
time-dependent  system.  The  implicit-explicit  method  was  rewritten  following  the 
methodology  of  [27]  with  the  hope  that  the  preconditioning  matrix  would  accelerate 
the  overall  convergence  of  the  scheme.  It  turned  out  that  this  idea  did  not  work 
well  for  the  I-E  scheme;  not  only  the  performance  did  not  improve  but  for  some 
runs  it  resulted  in  diverging  solutions.  This  behavior  may  be  attributed  to  two 
different  reasons:  first,  the  analysis  of  [27]  is  dependent  upon  the  finite-difference 
operator  and  the  results  quoted  in  [27]  do  not  carry  over  to  the  I-E  algorithm 
automatically,  and  second,  it  is  not  at  all  clear  that  the  idea  should  improve  time- 
dependent  solutions  in  the  first  place;  it  is  proven  to  work  well  for  steady-state 
solutions  but  there  does  not  exist  enough  evidence  that  the  acceleration  carries 
over  to  the  time-dependent  solutions. 

We  expect  that  there  are  many  possible  applications  of  the  I-E  approach  taken 
here.  In  some  of  them  (e.g.,  nearly  incompressible  flow)  an  obvious  competing  ap¬ 
proach  is  to  resolve  the  stiffness  problem  at  the  level  of  the  governing  equations  by 
deriving  a  new  system,  asymptotically  valid  in  the  appropriate  limit,  which  is  not 
stiff  and  developing  numerical  methods  for  the  limit  system  (e.g.,  the  incompress¬ 
ible  Euler  equations).  However,  it  should  be  noted  that  it  is  not  known  whether  or 
not  this  approach  is  valid  in  all  situations;  for  example,  in  the  case  of  magnetohy¬ 
drodynamics  in  the  limit  as  the  Alfven  number  approaches  zero,  it  does  not  appear 
to  be  possible  to  derive  a  limiting  set  of  equations,  see  [18].  Also,  it  is  usually  more 
natural  to  derive  well-posed  initial-boundary  value  problems  and  design  numerical 
boundary  conditions  for  the  original  set  of  equations  than  for  the  reduced  set  in 
the  appropriate  approximation.  Another  class  of  problems  for  which  our  approach 
is  an  obvious  candidate  arises  in  atmospheric  flows.  Here,  in  addition  to  the  Mach 
number  varying  widely,  low  amplitude  gravity  waves  are  present  and  restrict  the 
time  step  considerably;  eliminating  the  sound  waves  through  the  implementation  of 
reduced  equation  sets  can  lead  to  numerical  difficulties  in  specifying  open  boundary 
conditions. 
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The  application  of  shock  wave  -  boundary  layer  interaction,  which  we  were  not 
concerned  with  in  this  study,  requires  the  hybrid  scheme  to  transition  from  a  fully 
explicit  scheme  in  the  free  stream,  to  a  hybrid  scheme  at  the  top  of  the  boundary 
layer  and  possibly  to  a  fully  implicit  scheme  close  to  the  wall.  The  stability  results 
in  Part  II  indicate  that  the  transition  to  a  fully  implicit  scheme  in  a  two  dimensional 
flow  may  lead  to  instabilities.  These  instabilities  may  never  arise  in  this  application 
for  the  following  reasons:  the  flow  in  the  boundary  layer  near  the  wall  is  essentially 
one-dimensional  and  from  the  analysis  we  know  the  flow  has  to  be  fully  two  dimen¬ 
sional  for  the  instabilities  to  occur;  and  it  may  be  the  case  that  the  viscous  terms 
stabilize  the  scheme  in  this  region.  We  have  performed  preliminary  calculations 
using  a  viscous  extension  of  the  1-dimensional  method  from  Part  I  and  directional 
operator  splitting  for  shock  wave  -  boundary  layer  interactions;  the  results  indicate 
that  our  approach  will  be  useful  in  this  context.  The  implementation  of  a  viscous 
extension  of  the  unsplit  code  from  Part  II  will  be  a  topic  for  future  work. 

1.1.2  Chemically  Reacting  and  Dusty  Gas  Flows. 

The  objective  of  this  research  effort  is  the  efficient  and  accurate  computation  of 
reactive  multiphase  flows  at  both  high  and  low  Mach  numbers  (e.g.,  boundary 
layers,  combustion,  among  others).  We  believe  that  the  I— E  strategy,  when  coupled 
with  mesh  refinement  approaches  and  methods  for  the  stable  computation  of  stiff 
source  terms,  can  be  a  valuable  tool  in  this  field.  We  start  with  the  description  of 
preliminary  study  reported  in  [19]. 

The  purpose  of  the  study  in  [19]  was  to  further  advance  and  validate  the  I-E 
approach  by  investigating  the  effects  of  including  source  terms  representing  mul¬ 
tispecies  reacting  flows.  The  integration  scheme  for  the  sources  was  done  either 
fully  explicit  or  using  the  trapezoidal  rule,  but  in  both  cases  the  sources  were  not 
projected  onto  left  and  right  states  in  the  setup  of  initial  data  for  the  flux  function 
approximation  (or  Riemann  problem).  Several  ID  results  were  obtained;  converged 
solutions  for  variable  area  duct,  for  both  smooth  and  shocked  flows  are  demon¬ 
strated  using  an  extension  of  the  Newton  iteration  methodology  of  Collins  et  al 
[14].  Also,  an  appropriate  extension  of  the  unsplit  I-E  scheme  has  been  developed 
and  implemented  for  2D  reactive  flow  with  periodic  boundary  conditions;  the  new 
scheme  was  used  for  problems  analogous  to  the  nonreactive  results  discussed  above 
for  the  zero  Mach  number  limit. 

For  the  ID  calculations  we  find  that  the  performance  criteria  depend  on  the  degree 
of  source  term  stiffness.  In  particular,  a  fairly  low  CFL  maximum  is  required 
(that  is,  to  obtain  converged  solutions)  during  parts  of  the  calculations.  Unlike  the 
ID  calculations,  it  is  found  that  trapezoidal  rule  differencing  for  the  source  terms 
is  required  for  stability  for  the  2D  problem  solved.  In  view  of  these  results  we 
decided  to  explore  more  complex  treatments  of  stiff  source  terms  and,  especially, 
the  coupling  between  wave  speed  and  source  term  stiffness.  These  issues  were 
addressed  in  the  work  reported  in  [22]  and  described  next. 
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The  work  discussed  here  is  concerned  with  the  problem  of  source-term  stiffness 
which  arises  in,  e.g.,  reactive  flow  and/or  two-phase  flow.  We  studied  several  new 
innovative  approaches  due  to  Yee  et  al  [23]  [25]  [30],  Pember  [26]  and  the  present 
authors.  Several  ideas  were  then  applied  in  the  context  of  the  Godunov-like,  hybrid 
implicit-explicit  scheme  so  that  the  latter  methodology  was  successfully  extended 
to  treat  both  stiff  and  non-stiff,  nonequilibrium,  chemically  reacting  flow  fields. 
This  is  done  using  a  completely  unsplit  operator:  both  the  two-dimensional  and 
the  source  terms  integrations  axe  fully  coupled.  For  the  fully  explicit  mode  of  the 
scheme,  the  procedure  used  allows  for  the  formal  retention  of  second-order  accuracy 
of  the  numerical  scheme. 

The  stiff  and  nonstiff  solvers  developed  here  were  used  to  solve  two  unsteady  flow 
cases:  a  doubly  periodic,  unstable,  low  Mach  number  shear  layer  flow  and  the  sim¬ 
ulation  of  the  evolution  of  a  shock  wave  due  to  the  transverse  injection  of  fluid 
into  supersonic  flow.  In  both  flow  fields,  and  for  several  reaction  models  (i.e., 
frozen  air  flow,  air  with  only  vibrational  relaxation  of  molecules  and  a  complete 
chemically  reacting,  5-species  air  model),  the  numerical  results  prove  that  the  ex¬ 
tended  implicit-explicit  methodology  can  handle  complex  flow  problems  and  com¬ 
pute  high-resolution  details.  We  found  that  the  standard  non-stiff  approach  cannot 
handle  stiff  problems;  even  upon  starting  up  with  very  small  integration  steps,  the 
solution  becomes  unstable  and  eventually  breaks  down. 

We  have  also  shown  via  numerical  experimentation  that  after  the  resolution  of 
the  initial  large  gradients  in  stiff  flow  fields,  using  the  fully  explicit  mode  of  the 
solver,  the  solution  remains  stable  in  the  neighborhood  of  the  thermodynamical 
equilibrium  state  and  can  be  advanced  in  time  using  the  implicit-explicit  mode; 
although  the  I-E  mode  does  not  theoretically  provide  the  same  order  of  accuracy, 
it  may  be  used  due  to  its  relative  efficiency.  Therefore  the  stiff  I-E  solver  can  handle 
stiffness  arising  both  from  sources  and  wave  speeds. 

None  of  the  above  reports  contain  our  work  on  dust  modeling.  The  code  described 
here  contains  most  of  the  modules  needed  to  model  stiff  dusty  gas  flows.  Although 
similar  in  many  respects  to  the  problems  of  reactive  flow,  there  are  additional  nu¬ 
merical  difficulties  which  have  been  addressed  in  the  course  of  our  code  development 
work.  Future  tasks  are  to  optimize  the  code,  specifically  by  searching  for  and  im¬ 
plementing  faster  nonlinear  solvers  in  the  implicit  modes  of  the  predictor  and  the 
corrector  steps  (e.g.,  multigrid  methodologies,  see  [12]).  Also,  the  importance  of 
the  Riemann  problem  solver  (and  related  numerical  models)  needs  to  be  studied  in 
conjunction  with  the  solution  of  stiff,  close-to-equilibrium  thermodynamic  states. 
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1.2  NONCLASSICAL  GAS  DYNAMICS. 


1.2.1  Nonconvex  EOS. 

The  motivation  for  the  study  of  2D  unsteady  isentropic  gas  dynamics  (2x2  system) 
with  a  nonconvex  equation  of  state  (EOS),  as  well  as  the  results  of  the  work  that  was 
carried  out  during  1989-1991,  are  described  in  detail  in  [17],  [20].  A  brief  summary 
of  the  main  objectives  follows:  first,  to  simulate  the  a-e  phase  transition  in  Iron 
which  occurs  at  about  137  kilobars  and  has  been  extensively  studied;  second,  to 
test  several  variants  of  the  second— order  Godunov  scheme  in  the  simple  2x2  setting, 
but  for  materials  with  nonconvex  EOS. 

The  nature  of  the  above  mentioned  study  restricted  us  to  using  very  particular 
nonconvex  EOS’s.  We  used  two  models.  The  first  was  obtained  by  using  simple 
polytropic  gas  relations  such  as  p(v)  =  v 1 ,  where  p  is  the  pressure,  v  is  the  spe¬ 
cific  volume,  and  7  is  the  ratio  of  specific  heats,  and  joining  two  (locally)  convex 
disjoint  ’’pieces”  of  the  above  curve  by  adding  a  Cq°  function  so  that  the  resulting 
nonconvex  EOS  has  two  isolated  zeroes  in  the  second  derivative  and  the  function 
is  monotonically  decreasing.  The  second  model  was  based  on  experimental  data  of 
the  two  distinct  a—e  phases  of  Iron  along  the  isotherm  T  =  295°  K;  since  the  EOS 
curves  describing  the  two  phases  do  not  cross  in  the  transition  region,  we  carefully 
constructed  a  smooth  joining  (transition)  curve  which  consists  of  a  fifth  order  poly¬ 
nomial  with  second  order  contact  at  the  end  points.  As  a  consequence,  and  for 
both  models,  the  Riemann  problem  defined  for  any  two  initial  states  on  the  EOS 
curve  has  a  unique  solution,  composed  of  both  nonclassical  and  classical  waves, 
that  can  be  theoretically  predicted  (see,  e.g.,  [28])  and  numerically  constructed;  we 
point  out  that  such  solutions  are  complex,  tedious  to  program  and  very  expensive 
to  compute  exactly. 

The  results  of  the  above  mentioned  work  are  given  briefly  here.  First,  we  developed 
robust  and  accurate  numerical  schemes  for  the  solution  of  problems  described  by 
2x2  model  equations  with  a  nonconvex  EOS;  specifically,  we  found  that  the  gen¬ 
eral  methodology  of  Bell-Colella-Trangenstien  (BCT)  [4]  yielded  excellent  quality 
results  at  acceptable  expense  compared  to  other  Godunov-like  schemes.  Second, 
we  conducted  a  limited  study  of  self-similar,  oblique  shock  wave  reflections  in  two 
space  dimensions.  We  have  found  new  and  quite  exciting  phenomenology  resulting 
from  the  nonconvexity  that  induces  highly  nontrivial  2D  wave  interactions. 

The  work  summarized  above  has  been  extended  in  two  directions.  First,  we  ob¬ 
tained  more  results  in  the  2x2  context;  we  used  the  BCT  schemes  to  perform  new 
calculations  in  which  the  EOS  has  a  cusp.  The  latter  seems  to  more  realistically 
describe  phase  transitions  compared  to  the  smooth  curves  constructed  for  the  tran¬ 
sition  region  in  our  preceding  work  described  above.  The  results  are  extremely 
interesting,  especially  for  regular  to  Mach  transition,  and  they  clearly  indicate  the 
capturing  of  split  waves  as  predicted  by  the  theory  for  such  EOS  models.  More 
runs  and  analysis  are  required  for  the  better  understanding  of  this  phenomenon. 
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A  summary  of  this  ongoing  work  can  be  found  in  [21].  Second,  and  more  impor¬ 
tantly,  we  have  extended  the  BCT  algorithm  to  the  full  equations  of  gas  dynamics 
(3x3  system).  The  current  draft  of  E.  Erskine’s  Masters  thesis  [15]  contains  the 
fundamental  results  necessary  to  extend  the  BCT  algorithm  to  the  3x3  case  and  do 
realistic  material  modeling;  some  preliminary  numerical  results  in  1-space  dimen¬ 
sion  setting  were  presented  by  E.  Erskine  in  April  1993  at  a  meeting  in  Stonybrook, 
NY.  More  advanced  code  development  work  for  the  3x3  case  took  place  in  the  last 
few  months.  A  generic  version  of  a  BCT-like,  split,  two-spatial  dimensions  code  was 
written  and  results  for  a  convex  strong  shock  reflection  were  obtained  and  compared 
with  the  solutions  obtained  using  the  “usual”  second-order  Godunov  scheme;  the 
comparison  shows  excellent  agreement.  Currently,  several  nonconvex  EOS  are  being 
analytically  constructed  (a  nontrivial  issue  for  the  3x3  case  where  the  nonconvexity 
is  constructed  in  the  energy  equation  under  several  thermodynamical  constraints) 
and  they  will  be  tested  in  the  two-dimensional,  shock-reflection  setting  (some  very 
early  results  are  available).  Also,  we  are  surveying  the  literature  for  experimental 
data  that  could  assist  us  in  modeling  realistic  materials.  We  point  out  that  the 
results  for  the  nonconvex  3x3  case  are  new,  it  is  not  clear  how  the  physical  model 
translates  into  the  reflection  diagram  and  we  believe  that  a  large  number  of  runs 
and  exhaustive  analysis  will  be  required  to  verify  the  solutions  and  to  have  a  good 
feel  and  understanding  of  the  new  phenomenon. 

1.2.2  Multimaterial  Algorithms. 

For  a  description  of  the  main  issues  and  applications  for  multimaterial  extensions  of 
our  codes  see  [17].  During  the  contract  period,  it  became  apparent  that  our  method 
was  failing  for  the  case  of  a  strong  rarefaction  in  water.  Consequently,  considerable 
effort  has  gone  into  extending  our  analysis  and  fixing  the  code.  This  work  is  now 
complete  as  well  as  successfull  -  at  least  for  the  Tait  EOS.  Details  are  available  in 
the  current  draft  of  [9]. 
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1.3  INCOMPRESSIBLE  FLOWS  AND  THE  PROJECTION  METHOD. 


The  early  development  of  the  projection  method  is  discussed  in  [17].  Since  then,  J. 
Bell  and  colleagues  at  LLNL  worked  at  extending  the  code  to  treat  variable  density, 
3-D  flow  fields  and  also  combustion  in  the  zero  Mach  number  limit  [24].  The 
development  of  a  3-D  code  with  Adaptive  Mesh  Refinement  (AMR)  capabilities  is 
ongoing  at  LLNL.  After  completion,  it  will  be  possible  to  use  the  latter  to  simulate 
the  low  speed  part  of  fireball  rise.  In  what  follows  we  will  describe  some  preliminary 
studies  that  were  conducted  with  a  non-AMR  (i.e.,  uniform  mesh),  2-D  version  of 
the  code  that  was  provided  by  D.  Marcus  (LLNL). 

The  idea  is  to  numerically  model  the  cloud  rise  experiments  known  as  GEST  (Gas 
Explosive  Simulation  Technique,  see  [3]-[4])  by  first  solving  the  explosion  phase 
using  a  compressible  flow  solver  and,  after  some  appropriate  evolution  time,  as¬ 
suming  that  the  residual  flow  field  in  the  explosion  vicinity  is  incompressible,  carry 
out  the  solution  of  the  cloud  rise  phase  using  an  incompressible  flow  solver.  We 
can  proceed  as  follows.  A  spherical  explosion  was  modeled  incorporating  the  (ID) 
Godunov  compressible  code  and  using  the  input  data  produced  by  A.  Kuhl  (stoi¬ 
chiometric  Methane-Oxygen  mixture).  The  computation  was  carried  out  by  us  to 
the  physical  time  T  =  1.0  sec.,  for  which  the  leading  shock  wave  has  already  left  the 
computational  domain  (and  region  of  interest  for  the  incompressible  flow  model). 
We  remark  that  this  part  of  the  computation  was  done  on  a  relatively  fine  (and 
uniform)  mesh. 

The  output  data  obtained  by  the  latter  computation,  in  its  conservative  form,  is 
used  as  input  data  for  the  second  phase  of  the  computation.  First,  the  (spher¬ 
ical)  one- dimensional  conservative  variables  are  mapped  onto  a  two-dimensional, 
cylindrical  (r-z)  coordinate  system;  this  is  done  by  insisting  that  the  variables  re¬ 
main  conservative  with  respect  to  the  new  cell-volumes,  created  by  the  relatively 
coarse  cylindrical  mesh,  and  assigning  the  appropriate  (i.e.,  volumetric)  weight  to 
the  computed  variables  by  tracing  their  fine-mesh-volume  contributions  to  the  new 
cells  (in  which  they  are  fully  contained  or  partially  intersected).  The  incompressible 
flow  model  requires  only  the  density  and  velocity  fields.  However,  we  also  map  the 
composition  (i.e.,  product  gas/air  ratio)  and  simply  “convect”  this  variable  in  the 
evolving  incompressible  flow  field  in  order  to  extract  temperature  data  from  the 
local  density  and  composition  data. 

As  mentioned  above,  the  incompressible  solver  we  have  at  present  is  a  non-AMR,  2D 
version  of  the  second-order  projection  method.  The  computational  domain  chosen 
for  the  intended  preliminary  work  consists  of  256  by  512  zones,  modeling  a  physical 
domain  of  50m  (r)  by  100m  (z).  The  initial  data  consists  of  a  post-explosion  sphere 
of  radius  of  30m  mapped  onto  the  cylindrical  domain  and  centered  at  the  point  of 
ignition  of  the  explosive  device;  the  zones  outside  the  sphere  take  on  the  ambient 
air  conditions.  The  important  underlying  hypothesis  of  a  divergence  free  velocity 
field  for  the  incompressible  flow  model  is  (numerically)  imposed  on  the  initial  data 
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at  the  beginning  of  the  computation.  In  addition,  we  can  code  and  compute  the 
physical  locations  of  the  ten  probes  that  recorded  data  in  the  live  experiment,  tag 
the  (numerical)  cells  “containing”  these  probes,  and  record  time-dependent  data  at 
these  ten  locations.  Several  other  diagnostic  data  can  also  be  recorded  during  the 
’’cloud  rise”  computational  phase.  A  complete  run,  one  which  advances  the  flow 
from  T  =  1.0  sec.  (end  of  compressible  flow  solver)  to  about  T  =  4.0  sec.  (the 
data  recordings  of  the  successful  experiment  terminated  at  about  T  =  3.8  sec.),  is 
expected  to  take  about  800  time  steps  at  a  cost  of  about  one  CPU  hour  on  a  Cray 
YMP. 

The  results  would  consist  of  the  time  dependent  density,  velocity,  vorticity  and 
composition  fields.  The  vorticity  field  will  capture  the  roll  up  of  vortices  and  their 
motion  upwards,  based  on  previous  work  with  this  code  technology.  A  full  analysis 
of  such  a  calculation  must  await  the  results.  We  are  now  completely  prepared  to 
proceed  with  this  effort. 
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SECTION  2 
CONCLUSIONS 


(1)  The  implicit-explicit  methodology  has  been  extended  to  near  zero  Mach 
number  situations  and  accurate,  resolved  results  have  been  computed. 
Source  term  capability  -  stiff  and  nonstiff  -  has  been  successfully  imple¬ 
mented  in  the  context  of  realistic  physics  modeling. 

(2)  The  variable  density  projection  algorithm  is  a  highly  robust,  accurate  and 
efficient  method  for  problems  which  axe  essentially  zero  Mach  number  at 
all  times  of  interest,  e.g.,  cloud  rise.  It  is  the  method  of  choice  for  such 
situations. 

(3)  The  multimaterial  algorithm  is  now  available  for  water  using  the  Tait  EOS. 
More  generally,  the  algorithm  is  now  being  used  by  many  groups  with  great 
success. 

(4)  It  has  been  shown  that  materials  characterized  by  nonclassical  EOS  exhibit 
solution  phenomenology  which  can  be  profoundly  different  from  convex 
EOS  gases.  The  implications  of  this  are  quite  profound  for  the  numerical 
prediction  of  shock  and  blast  wave  loading;  it  is  necessary  that  the  under¬ 
lying  code  integrator  be  capable  of  capturing  all  of  the  physically  relevant 
waves  and  do  so  if  and  only  if  they  are  actually  present  in  the  mathematical 
solution. 
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SECTION  3 

RECOMMENDATIONS 

(1)  The  implicit-explicit  code  remains  one  of  two  possibilities  for  computing 
problems  with  inherent  wave  speed  stiffness,  along  with  hybrids  of  pro¬ 
jection  codes  and  explicit  upwind  schemes.  Further  research  is  essential 
if  endusers  are  to  have  confidence  in  the  results  from  simulations  of  prob¬ 
lems  such  as  fireball  rise  where  the  Mach  number  can  increase  substantially 
during  the  course  of  the  calculation. 

(2)  The  explicit  codes  for  gas  dynamics  are  now  a  fairly  mature  technology.  Our 
research  on  multimaterial  extensions  and  more  complicated  conservation 
laws  leads  to  the  conclusion  that  code  development  efforts  in  this  direction 
would  payoff  in  greatly  improved  predictive  capabilities  for  applications 
such  as  blast  wave  loading  of  soil,  reinforced  targets,  etc.  and  advanced 
HE  weapon  design  (including  interior  ballistics). 

(3)  The  zero  Mach  number  projection  methodology  is  advancing  rapidly  and 
it  appears  that  it’s  extension  to  long-time  transport  at  large  lenghth  scales 
(including  mesoscale)  will  pose  no  theoretical  difficulties.  In  view  of  the 
relatively  low  level  of  sophistication  in  the  numerical  modelling  in  currently 
available  codes  for  this  problem  domain,  a  large  effort  here  promises  even 
larger  payoffs. 
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